COVID-19 US

Prepare county level data

Read and format personality data

df_us_pers <- read_csv('timeseries_usa_county_march1_april_09.csv')

df_us_pers <- df_us_pers %>% select(countyfips, open, sci, extra, agree, stabil) %>% 
  mutate(stabil = 6-stabil) %>%
  dplyr::rename(county_fips = countyfips,
         pers_o = open, 
         pers_c = sci,
         pers_e = extra,
         pers_a = agree,
         pers_n = stabil) %>% 
  distinct() %>%
  mutate(county_fips = as.character(county_fips))

df_us_pers

Read and format prevalence data

df_us_prev <- read_csv('USA_timeseries_prep_2005.csv')

df_us_prev <- df_us_prev %>% 
  select(fips, date, rate) %>% 
  mutate(date = as.Date(date, "%d%b%Y")) %>% 
  rename(county_fips = fips, 
         rate_day = rate) %>%
  mutate(county_fips = as.character(county_fips))

df_us_prev

NA

Read and format county level controls

df_us_ctrl <- read.csv('controls_US.csv')

df_us_ctrl <- df_us_ctrl %>% select(-county_name) %>% 
  rename(county_fips = county) %>%
  mutate(county_fips = as.character(county_fips))

df_us_ctrl

NA

Read and format social distancing data FB

fb_files <- list.files('../FB Data/US individual files/Mobility/',
                       '*.csv', full.names = T)

df_us_socdist <- fb_files %>% 
  map(read_csv) %>% bind_rows()

df_us_socdist <- df_us_socdist %>%
  select(-age_bracket, -gender, -baseline_name, -baseline_type, -polygon_name) %>%
  rename(date = ds,
         county_fips = polygon_id,
         socdist_tiles = all_day_bing_tiles_visited_relative_change,
         socdist_single_tile = all_day_ratio_single_tile_users) %>%
  mutate(county_fips = as.character(county_fips))

df_us_socdist

NA

Merge data

# create sequence of dates
date_sequence <- seq.Date(min(df_us_prev$date),
                          max(df_us_prev$date), 1)
                     
# create data frame with time sequence
df_dates = data.frame(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')

# join data frames 
df_us_prev <- df_us_prev %>%
  plyr::join(df_us_ctrl, by='county_fips') %>% 
  plyr::join(df_us_pers, by='county_fips') %>%
  merge(df_dates, by='date') %>% 
  arrange(county_fips, date)

df_us_prev

# create sequence of dates
date_sequence <- seq.Date(min(df_us_socdist$date),
                          max(df_us_socdist$date), 1)
                     
# create data frame with time sequence
df_dates = data.frame(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')

# join data frames 
df_us_socdist <- df_us_socdist %>%
  plyr::join(df_us_ctrl, by='county_fips') %>% 
  plyr::join(df_us_pers, by='county_fips') %>%
  merge(df_dates, by='date') %>% 
  arrange(county_fips, date)

fips_complete <- df_us_socdist %>% 
  group_by(county_fips) %>% 
  summarize(n = n()) %>% 
  filter(! n<max(n)) %>% .$county_fips

df_us_socdist <- df_us_socdist %>%
  filter(county_fips %in% fips_complete)

df_us_socdist

Control for weekend effect

df_us_loess <- df_us_socdist %>% 
  mutate(weekday = format(date, '%u')) %>% 
  filter(!weekday %in% c('6','7')) %>% 
  split(.$county_fips) %>%
  map(~ loess(socdist_single_tile ~ time, data = .)) %>%
  map(predict, 1:max(df_us_socdist$time)) %>% 
  bind_rows() %>% 
  gather(key = 'county_fips', value = 'loess') %>% 
  group_by(county_fips) %>% 
  mutate(time = row_number())

df_us_socdist <- df_us_socdist %>% merge(df_us_loess, by=c('county_fips', 'time')) %>% 
  mutate(weekday = format(date, '%u')) %>% 
  mutate(socdist_single_tile_clean = ifelse(weekday %in% c('6','7'), loess,
                                            socdist_single_tile)) %>%
  arrange(county_fips, time) %>% 
  select(-weekday)

df_us_socdist <- df_us_socdist %>% drop_na() %>% mutate(time = time-1)

Plot prevalence over time

df_us_prev %>% sample_n(20000) %>%
  ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall prevalence over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_us_prev %>% 
  mutate(prev_tail = cut(.[[i]], 
                         breaks = c(-Inf, quantile(.[[i]], 0.05, na.rm=T), 
                                    quantile(.[[i]], 0.95, na.rm=T), Inf),
                         labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(prev_tail != 'center') %>%
  ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~prev_tail) + 
  theme(legend.position="none") + 
  ggtitle(i)

print(gg)
}

Plot social distancing single tile visited

df_us_socdist %>% sample_n(10000) %>%
  ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing (single tile) over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_us_socdist %>% 
  mutate(dist_tail = cut(.[[i]], 
                         breaks = c(-Inf, quantile(.[[i]], 0.05, na.rm = T), 
                                    quantile(.[[i]], 0.95, na.rm = T), Inf),
                         labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(dist_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~dist_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

df_us_socdist <- df_us_socdist %>% mutate(socdist_single_tile = socdist_single_tile_clean) %>% 
  select(-loess, -socdist_single_tile_clean)

Correlations

df_us_prev %>% select(-time, -date) %>% 
  group_by(county_fips) %>%
  summarize_if(is.numeric, mean) %>% 
  select(-county_fips) %>%
  cor(use='pairwise.complete.obs') %>% 
  round(3) %>% as.data.frame()

df_us_socdist %>% select(-time, -date) %>% 
  group_by(county_fips) %>%
  summarize_if(is.numeric, mean) %>% 
  select(-county_fips) %>%
  cor(use='pairwise.complete.obs') %>% 
  round(3) %>% as.data.frame()

NA

Rescale Data

lvl2_scaled <- df_us_prev %>% 
  select(-time, -date, -rate_day) %>% 
  distinct() %>% 
  mutate_at(vars(-county_fips), scale)

lvl1_scaled <- df_us_prev %>% select(county_fips, time, rate_day)

df_us_prev_scaled <- plyr::join(lvl1_scaled, lvl2_scaled, by = 'county_fips') 

lvl2_scaled <- df_us_socdist %>% 
  select(-time, -date, -socdist_tiles, -socdist_single_tile) %>% 
  distinct() %>% 
  mutate_at(vars(-county_fips), scale)

lvl1_scaled <- df_us_socdist %>% select(county_fips, time, socdist_single_tile)

df_us_socdist_scaled <- plyr::join(lvl1_scaled, lvl2_scaled, by = 'county_fips') 

Predict Prevalence

Extract first day of covid outbreak

# get onset day
df_us_onset_prev <- df_us_prev_scaled %>% 
  group_by(county_fips) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  summarize(onset_prev = min(time)) %>%
  mutate(county_fips = as.character(county_fips))
  
# merge with county data
df_us_onset_prev <- df_us_prev_scaled %>% 
  select(-time, -rate_day) %>%
  distinct() %>% 
  mutate(county_fips = as.character(county_fips)) %>%
  left_join(df_us_onset_prev, by = 'county_fips')

# handle censored data
df_us_onset_prev <- df_us_onset_prev %>% 
  mutate(event = ifelse(is.na(onset_prev), 0, 1)) %>% 
  mutate(onset_prev = replace_na(onset_prev, as.numeric(diff(range(df_us_prev$date)))+1))

Extract slopes

# cut time series before onset
df_us_prev_scaled <- df_us_prev_scaled %>% 
  group_by(county_fips) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  mutate(time = time-min(time)+1) %>%
  ungroup() %>%
  filter(time <= 30) %>%
  select(-rate_cs)

# drop counties with little data
df_us_prev_scaled <- df_us_prev_scaled %>%
  group_by(county_fips) %>%
  filter(n() == 30) %>%
  ungroup()

# extract slope prevalence
df_us_slope_prev <- df_us_prev_scaled %>% split(.$county_fips) %>% 
  map(~ lm(rate_day ~ time, data = .)) %>%
  map(coef) %>% 
  map_dbl('time') %>% 
  as.data.frame() %>% 
  rownames_to_column('county_fips') %>% 
  rename(slope_prev = '.')

# merge with county data
df_us_slope_prev <- df_us_prev_scaled %>% 
  select(-time, -rate_day) %>%
  distinct() %>% 
  inner_join(df_us_slope_prev, by = 'county_fips') %>%
  drop_na()

Explore distributions

df_us_onset_prev %>% ggplot(aes(onset_prev)) + geom_histogram(bins = 100)

df_us_slope_prev %>% ggplot(aes(slope_prev)) + geom_histogram(bins = 100)

Predict COVID onset with time-to-event regression

# predict onset from personality
cox_onset_prev <- coxph(Surv(onset_prev, event) ~ 
                          pers_o + pers_c + pers_e + pers_a + pers_n, 
                        data = df_us_onset_prev)
cox_onset_prev %>% summary()

Call:
coxph(formula = Surv(onset_prev, event) ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n, data = df_us_onset_prev)

  n= 2469, number of events= 2429 
   (440 observations deleted due to missingness)

           coef exp(coef) se(coef)      z Pr(>|z|)    
pers_o  0.25801   1.29435  0.01953 13.214  < 2e-16 ***
pers_c  0.01242   1.01250  0.02440  0.509   0.6107    
pers_e  0.04044   1.04127  0.01961  2.062   0.0392 *  
pers_a  0.02804   1.02844  0.02404  1.167   0.2434    
pers_n -0.13832   0.87082  0.02214 -6.248 4.17e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    1.2944     0.7726    1.2458    1.3448
pers_c    1.0125     0.9877    0.9652    1.0621
pers_e    1.0413     0.9604    1.0020    1.0821
pers_a    1.0284     0.9723    0.9811    1.0781
pers_n    0.8708     1.1483    0.8338    0.9094

Concordance= 0.658  (se = 0.006 )
Likelihood ratio test= 280.5  on 5 df,   p=<2e-16
Wald test            = 282.1  on 5 df,   p=<2e-16
Score (logrank) test = 280.1  on 5 df,   p=<2e-16

# predict onset from personality with controls
cox_onset_prev_ctrl <- coxph(Surv(onset_prev, event) ~ 
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_distance + republican + medage + male + popdens +
                               manufact + tourism + academics + medinc + physician_pc,
                             data = df_us_onset_prev)
cox_onset_prev_ctrl %>% summary()

Call:
coxph(formula = Surv(onset_prev, event) ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n + airport_distance + republican + medage + 
    male + popdens + manufact + tourism + academics + medinc + 
    physician_pc, data = df_us_onset_prev)

  n= 2464, number of events= 2424 
   (445 observations deleted due to missingness)

                      coef exp(coef)  se(coef)       z Pr(>|z|)    
pers_o            0.219404  1.245334  0.024013   9.137  < 2e-16 ***
pers_c            0.074751  1.077615  0.026097   2.864 0.004179 ** 
pers_e            0.016376  1.016510  0.020873   0.785 0.432729    
pers_a            0.042856  1.043787  0.026698   1.605 0.108444    
pers_n            0.001372  1.001373  0.025862   0.053 0.957701    
airport_distance -0.293113  0.745938  0.029875  -9.811  < 2e-16 ***
republican       -0.181505  0.834014  0.025588  -7.093 1.31e-12 ***
medage           -0.277857  0.757405  0.022256 -12.484  < 2e-16 ***
male             -0.145115  0.864923  0.028091  -5.166 2.39e-07 ***
popdens          -0.090498  0.913477  0.024988  -3.622 0.000293 ***
manufact          0.046112  1.047192  0.024014   1.920 0.054827 .  
tourism           0.068927  1.071358  0.021260   3.242 0.001186 ** 
academics         0.052582  1.053989  0.036863   1.426 0.153748    
medinc            0.308101  1.360838  0.032191   9.571  < 2e-16 ***
physician_pc     -0.107555  0.898027  0.022305  -4.822 1.42e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                 exp(coef) exp(-coef) lower .95 upper .95
pers_o              1.2453     0.8030    1.1881    1.3053
pers_c              1.0776     0.9280    1.0239    1.1342
pers_e              1.0165     0.9838    0.9758    1.0590
pers_a              1.0438     0.9580    0.9906    1.0999
pers_n              1.0014     0.9986    0.9519    1.0534
airport_distance    0.7459     1.3406    0.7035    0.7909
republican          0.8340     1.1990    0.7932    0.8769
medage              0.7574     1.3203    0.7251    0.7912
male                0.8649     1.1562    0.8186    0.9139
popdens             0.9135     1.0947    0.8698    0.9593
manufact            1.0472     0.9549    0.9990    1.0977
tourism             1.0714     0.9334    1.0276    1.1169
academics           1.0540     0.9488    0.9805    1.1330
medinc              1.3608     0.7348    1.2776    1.4495
physician_pc        0.8980     1.1136    0.8596    0.9382

Concordance= 0.735  (se = 0.005 )
Likelihood ratio test= 998.1  on 15 df,   p=<2e-16
Wald test            = 1107  on 15 df,   p=<2e-16
Score (logrank) test = 1225  on 15 df,   p=<2e-16

Predict prevalence slopes with linear models

# predict slopes from personality
lm_slope_prev <- lm(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_us_slope_prev)
lm_slope_prev %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, data = df_us_slope_prev)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.08820 -0.03714 -0.02402 -0.00237  2.96494 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.048019   0.002665  18.016  < 2e-16 ***
pers_o      -0.006137   0.003041  -2.018  0.04371 *  
pers_c      -0.002526   0.003685  -0.685  0.49311    
pers_e      -0.002132   0.003058  -0.697  0.48581    
pers_a       0.010416   0.003754   2.774  0.00558 ** 
pers_n      -0.005947   0.003448  -1.725  0.08468 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1265 on 2254 degrees of freedom
Multiple R-squared:  0.01033,   Adjusted R-squared:  0.008138 
F-statistic: 4.707 on 5 and 2254 DF,  p-value: 0.0002796

# predict slopes from personality with controls
lm_slope_prev_ctrl <- lm(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n + 
                           airport_distance + republican + medage + male + popdens + 
                           manufact + tourism + academics + medinc + physician_pc,
                         data = df_us_slope_prev)
lm_slope_prev_ctrl %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + airport_distance + republican + medage + male + 
    popdens + manufact + tourism + academics + medinc + physician_pc, 
    data = df_us_slope_prev)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.22382 -0.03886 -0.01734  0.00683  2.79539 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       0.046634   0.002760  16.894  < 2e-16 ***
pers_o           -0.001848   0.003533  -0.523  0.60097    
pers_c           -0.002344   0.003688  -0.636  0.52502    
pers_e           -0.001526   0.003058  -0.499  0.61775    
pers_a            0.010595   0.003873   2.735  0.00628 ** 
pers_n           -0.003037   0.003719  -0.817  0.41420    
airport_distance  0.005851   0.003655   1.601  0.10952    
republican       -0.013491   0.003300  -4.088 4.50e-05 ***
medage           -0.008880   0.002915  -3.047  0.00234 ** 
male              0.008961   0.003493   2.566  0.01036 *  
popdens           0.007500   0.002733   2.744  0.00611 ** 
manufact          0.010568   0.003207   3.295  0.00100 ** 
tourism           0.001337   0.003518   0.380  0.70401    
academics        -0.011711   0.005040  -2.323  0.02025 *  
medinc            0.012500   0.004258   2.936  0.00336 ** 
physician_pc      0.013890   0.002918   4.760 2.06e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1242 on 2244 degrees of freedom
Multiple R-squared:  0.05072,   Adjusted R-squared:  0.04437 
F-statistic: 7.993 on 15 and 2244 DF,  p-value: < 2.2e-16

CRF predicting slopes

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_slope_prev <- cforest(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n + 
                           airport_distance + republican + medage + male + popdens + 
                           manufact + tourism + academics + medinc + physician_pc,
                           data = df_us_slope_prev, 
                         controls = ctrls)

crf_slope_prev_varimp <- varimp(crf_slope_prev, nperm = 1)
crf_slope_prev_varimp_cond <- varimp(crf_slope_prev, conditional = T, nperm = 1)

crf_slope_prev_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_slope_prev_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

Predict Social Distancing

Change point analysis

# keep only counties with full data
fips_complete <- df_us_socdist_scaled %>% 
  group_by(county_fips) %>% 
  summarize(n = n()) %>% 
  filter(n==max(.$n)) %>% 
  .$county_fips

# run changepoint analysis
df_us_socdist_cpt_results <- df_us_socdist_scaled %>% 
  select(county_fips, socdist_single_tile) %>%
  filter(county_fips %in% fips_complete) %>% 
  split(.$county_fips) %>%
  map(~ cpt.meanvar(as.vector(.$socdist_single_tile),
                    #penalty = 'Asymptotic',
                    class=TRUE,
                    param.estimates=TRUE,
                    Q=1,
                    test.stat = 'Normal'))

# calculate change point
df_us_socdist_cpt_day <- df_us_socdist_cpt_results %>% 
  map(cpts) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(cpt_day_socdist = '.') %>%
  rownames_to_column('county_fips')

# calculate mean differences
df_us_socdist_cpt_mean_diff <- df_us_socdist_cpt_results %>% 
  map(param.est) %>% 
  map(~ .$mean) %>% 
  map(~ .[2]-.[1]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(mean_diff_socdist = '.') %>%
  rownames_to_column('county_fips')

# calculate varaince differences
df_us_socdist_cpt_var_diff <- df_us_socdist_cpt_results %>% 
  map(param.est) %>% 
  map(~ .$variance) %>% 
  map(~ .[2]-.[1]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(var_diff_socdist = '.') %>%
  rownames_to_column('county_fips')

# merge with county data
df_us_cpt_socdist <- df_us_socdist_scaled %>% 
  select(-time, -socdist_single_tile) %>%
  distinct() %>% 
  left_join(df_us_socdist_cpt_day, by='county_fips') %>%
  left_join(df_us_socdist_cpt_mean_diff, by='county_fips') %>%
  left_join(df_us_socdist_cpt_var_diff, by='county_fips') %>%
  left_join(select(df_us_onset_prev, county_fips, onset_prev), by='county_fips') %>%
  left_join(select(df_us_slope_prev, county_fips, slope_prev), by='county_fips') 

# handle censored data
df_us_cpt_socdist <- df_us_cpt_socdist %>% 
  mutate(cpt_day_socdist = ifelse(is.na(cpt_day_socdist), as.numeric(diff(range(df_us$date))), cpt_day_socdist)) %>% 
  mutate(event = ifelse(cpt_day_socdist >= 60, 0, 1))

df_us_cpt_socdist$cpt_day_socdist %>% hist()

df_us_cpt_socdist$mean_diff_socdist %>% hist()

df_us_cpt_socdist$var_diff_socdist %>% hist()

for(i in head(df_us_socdist_cpt_results, 5)){
  plot(i)
}

NA

Predicting change points with time-to-event regression

# predict hazard from personality
cox_cpt_socdist <- coxph(Surv(cpt_day_socdist, event) ~ 
                           pers_o + pers_c + pers_e + pers_a + pers_n, 
                  data = df_us_cpt_socdist)
cox_cpt_socdist %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n, data = df_us_cpt_socdist)

  n= 2380, number of events= 2379 

           coef exp(coef) se(coef)      z Pr(>|z|)    
pers_o -0.15153   0.85939  0.02244 -6.754 1.44e-11 ***
pers_c -0.02846   0.97194  0.02681 -1.062  0.28843    
pers_e  0.02941   1.02985  0.02358  1.247  0.21226    
pers_a -0.08052   0.92264  0.02767 -2.910  0.00361 ** 
pers_n  0.08166   1.08508  0.02678  3.050  0.00229 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    0.8594     1.1636    0.8224    0.8980
pers_c    0.9719     1.0289    0.9222    1.0244
pers_e    1.0299     0.9710    0.9833    1.0786
pers_a    0.9226     1.0838    0.8739    0.9741
pers_n    1.0851     0.9216    1.0296    1.1435

Concordance= 0.582  (se = 0.008 )
Likelihood ratio test= 107.2  on 5 df,   p=<2e-16
Wald test            = 109.7  on 5 df,   p=<2e-16
Score (logrank) test = 109.3  on 5 df,   p=<2e-16

# predict hazard from personality with controls
cox_cpt_socdist_ctrl <- coxph(Surv(cpt_day_socdist, event) ~ 
                                pers_o + pers_c + pers_e + pers_a + pers_n + 
                           airport_distance + republican + medage + male + popdens + 
                           manufact + tourism + academics + medinc + physician_pc,
                  data = df_us_cpt_socdist)
cox_cpt_socdist_ctrl %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n + airport_distance + republican + 
    medage + male + popdens + manufact + tourism + academics + 
    medinc + physician_pc, data = df_us_cpt_socdist)

  n= 2380, number of events= 2379 

                      coef exp(coef)  se(coef)      z Pr(>|z|)    
pers_o           -0.094909  0.909456  0.026471 -3.585 0.000337 ***
pers_c           -0.024506  0.975792  0.027403 -0.894 0.371167    
pers_e            0.035056  1.035678  0.023966  1.463 0.143539    
pers_a           -0.091709  0.912370  0.029044 -3.158 0.001591 ** 
pers_n            0.071328  1.073933  0.029010  2.459 0.013942 *  
airport_distance  0.049305  1.050541  0.024199  2.037 0.041606 *  
republican        0.046018  1.047093  0.025565  1.800 0.071850 .  
medage            0.039634  1.040430  0.020516  1.932 0.053374 .  
male              0.008525  1.008561  0.021596  0.395 0.693038    
popdens          -0.130542  0.877620  0.038711 -3.372 0.000746 ***
manufact          0.102368  1.107792  0.023238  4.405 1.06e-05 ***
tourism          -0.073016  0.929586  0.024515 -2.978 0.002898 ** 
academics         0.167903  1.182822  0.040263  4.170 3.04e-05 ***
medinc           -0.086244  0.917370  0.033798 -2.552 0.010718 *  
physician_pc     -0.020033  0.980166  0.023267 -0.861 0.389235    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                 exp(coef) exp(-coef) lower .95 upper .95
pers_o              0.9095     1.0996    0.8635    0.9579
pers_c              0.9758     1.0248    0.9248    1.0296
pers_e              1.0357     0.9656    0.9882    1.0855
pers_a              0.9124     1.0960    0.8619    0.9658
pers_n              1.0739     0.9312    1.0146    1.1368
airport_distance    1.0505     0.9519    1.0019    1.1016
republican          1.0471     0.9550    0.9959    1.1009
medage              1.0404     0.9611    0.9994    1.0831
male                1.0086     0.9915    0.9668    1.0522
popdens             0.8776     1.1394    0.8135    0.9468
manufact            1.1078     0.9027    1.0585    1.1594
tourism             0.9296     1.0757    0.8860    0.9753
academics           1.1828     0.8454    1.0931    1.2799
medinc              0.9174     1.0901    0.8586    0.9802
physician_pc        0.9802     1.0202    0.9365    1.0259

Concordance= 0.6  (se = 0.008 )
Likelihood ratio test= 184.8  on 15 df,   p=<2e-16
Wald test            = 166.6  on 15 df,   p=<2e-16
Score (logrank) test = 171.9  on 15 df,   p=<2e-16

# predict hazard from personality with controls
cox_cpt_socdist_ctrl2 <- coxph(Surv(cpt_day_socdist, event) ~ 
                                 pers_o + pers_c + pers_e + pers_a + pers_n + 
                           airport_distance + republican + medage + male + popdens + 
                           manufact + tourism + academics + medinc + physician_pc +
                           onset_prev + slope_prev ,
                  data = df_us_cpt_socdist)
cox_cpt_socdist_ctrl2 %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n + airport_distance + republican + 
    medage + male + popdens + manufact + tourism + academics + 
    medinc + physician_pc + onset_prev + slope_prev, data = df_us_cpt_socdist)

  n= 2196, number of events= 2195 
   (184 observations deleted due to missingness)

                      coef exp(coef)  se(coef)      z Pr(>|z|)    
pers_o           -0.074555  0.928157  0.028688 -2.599 0.009354 ** 
pers_c           -0.025533  0.974790  0.029119 -0.877 0.380560    
pers_e            0.036051  1.036709  0.026042  1.384 0.166254    
pers_a           -0.079471  0.923605  0.031369 -2.533 0.011295 *  
pers_n            0.081537  1.084953  0.031671  2.574 0.010040 *  
airport_distance  0.051080  1.052407  0.026626  1.918 0.055057 .  
republican        0.013368  1.013457  0.026805  0.499 0.617996    
medage            0.007390  1.007418  0.022039  0.335 0.737375    
male              0.001792  1.001793  0.023217  0.077 0.938489    
popdens          -0.293911  0.745343  0.068400 -4.297 1.73e-05 ***
manufact          0.118303  1.125585  0.024685  4.792 1.65e-06 ***
tourism          -0.073760  0.928894  0.026138 -2.822 0.004773 ** 
academics         0.156968  1.169958  0.042147  3.724 0.000196 ***
medinc           -0.037278  0.963408  0.035395 -1.053 0.292256    
physician_pc     -0.020296  0.979909  0.024747 -0.820 0.412136    
onset_prev        0.003293  1.003299  0.003139  1.049 0.294063    
slope_prev       -0.869099  0.419329  0.239187 -3.634 0.000280 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                 exp(coef) exp(-coef) lower .95 upper .95
pers_o              0.9282     1.0774    0.8774    0.9818
pers_c              0.9748     1.0259    0.9207    1.0320
pers_e              1.0367     0.9646    0.9851    1.0910
pers_a              0.9236     1.0827    0.8685    0.9822
pers_n              1.0850     0.9217    1.0197    1.1544
airport_distance    1.0524     0.9502    0.9989    1.1088
republican          1.0135     0.9867    0.9616    1.0681
medage              1.0074     0.9926    0.9648    1.0519
male                1.0018     0.9982    0.9572    1.0484
popdens             0.7453     1.3417    0.6518    0.8523
manufact            1.1256     0.8884    1.0724    1.1814
tourism             0.9289     1.0765    0.8825    0.9777
academics           1.1700     0.8547    1.0772    1.2707
medinc              0.9634     1.0380    0.8988    1.0326
physician_pc        0.9799     1.0205    0.9335    1.0286
onset_prev          1.0033     0.9967    0.9971    1.0095
slope_prev          0.4193     2.3848    0.2624    0.6701

Concordance= 0.605  (se = 0.008 )
Likelihood ratio test= 189.8  on 17 df,   p=<2e-16
Wald test            = 172.7  on 17 df,   p=<2e-16
Score (logrank) test = 175  on 17 df,   p=<2e-16

Linear models predicting mean differences

lm_meandiff_socdist <- lm(mean_diff_socdist ~ 
                            pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_us_cpt_socdist)
lm_meandiff_socdist %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, data = df_us_cpt_socdist)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.083345 -0.018086 -0.002951  0.014510  0.175879 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.0921830  0.0005616 164.145  < 2e-16 ***
pers_o       0.0106893  0.0006089  17.554  < 2e-16 ***
pers_c      -0.0048685  0.0007619  -6.390 1.99e-10 ***
pers_e       0.0042813  0.0006264   6.835 1.04e-11 ***
pers_a      -0.0042184  0.0007778  -5.424 6.43e-08 ***
pers_n      -0.0030490  0.0007107  -4.290 1.86e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0274 on 2374 degrees of freedom
Multiple R-squared:  0.2192,    Adjusted R-squared:  0.2176 
F-statistic: 133.3 on 5 and 2374 DF,  p-value: < 2.2e-16

lm_meandiff_socdist_ctrl <- lm(mean_diff_socdist ~ 
                                 pers_o + pers_c + pers_e + pers_a + pers_n + 
                              airport_distance + republican + medage + male + popdens + 
                              manufact + tourism + academics + medinc + physician_pc,
                            data = df_us_cpt_socdist)
lm_meandiff_socdist_ctrl %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + airport_distance + republican + medage + male + 
    popdens + manufact + tourism + academics + medinc + physician_pc, 
    data = df_us_cpt_socdist)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.108336 -0.011582 -0.001643  0.010082  0.098895 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)       0.0921830  0.0003778 243.977  < 2e-16 ***
pers_o            0.0020818  0.0004813   4.325 1.59e-05 ***
pers_c           -0.0014328  0.0005213  -2.749 0.006028 ** 
pers_e            0.0015901  0.0004285   3.711 0.000211 ***
pers_a           -0.0008177  0.0005497  -1.488 0.137004    
pers_n            0.0044056  0.0005238   8.411  < 2e-16 ***
airport_distance -0.0023949  0.0004084  -5.865 5.13e-09 ***
republican       -0.0092408  0.0004570 -20.222  < 2e-16 ***
medage            0.0020335  0.0004009   5.072 4.23e-07 ***
male              0.0004270  0.0004170   1.024 0.305979    
popdens           0.0046320  0.0004051  11.436  < 2e-16 ***
manufact          0.0004644  0.0004546   1.021 0.307147    
tourism           0.0033340  0.0004515   7.385 2.11e-13 ***
academics         0.0051900  0.0007585   6.842 9.91e-12 ***
medinc            0.0132681  0.0006282  21.119  < 2e-16 ***
physician_pc     -0.0001920  0.0004176  -0.460 0.645750    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01843 on 2364 degrees of freedom
Multiple R-squared:  0.6481,    Adjusted R-squared:  0.6458 
F-statistic: 290.2 on 15 and 2364 DF,  p-value: < 2.2e-16

lm_meandiff_socdist_ctrl2 <- lm(mean_diff_socdist ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n + 
                              airport_distance + republican + medage + male + popdens + 
                              manufact + tourism + academics + medinc + physician_pc +
                              onset_prev + slope_prev ,
                            data = df_us_cpt_socdist)
lm_meandiff_socdist_ctrl2 %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + airport_distance + republican + medage + male + 
    popdens + manufact + tourism + academics + medinc + physician_pc + 
    onset_prev + slope_prev, data = df_us_cpt_socdist)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.111040 -0.011493 -0.001667  0.010145  0.075483 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)       1.194e-01  3.741e-03  31.914  < 2e-16 ***
pers_o            1.641e-03  5.132e-04   3.198 0.001405 ** 
pers_c           -1.520e-03  5.382e-04  -2.824 0.004782 ** 
pers_e            2.020e-03  4.492e-04   4.498 7.23e-06 ***
pers_a           -1.349e-03  5.762e-04  -2.342 0.019271 *  
pers_n            4.468e-03  5.479e-04   8.155 5.83e-16 ***
airport_distance -1.445e-03  4.359e-04  -3.316 0.000928 ***
republican       -8.186e-03  4.740e-04 -17.271  < 2e-16 ***
medage            3.138e-03  4.193e-04   7.485 1.03e-13 ***
male              1.172e-03  4.355e-04   2.691 0.007180 ** 
popdens           7.177e-03  9.791e-04   7.331 3.21e-13 ***
manufact          3.451e-04  4.655e-04   0.741 0.458640    
tourism           2.675e-03  4.683e-04   5.711 1.28e-08 ***
academics         4.908e-03  7.672e-04   6.397 1.93e-10 ***
medinc            1.188e-02  6.481e-04  18.333  < 2e-16 ***
physician_pc     -2.821e-04  4.264e-04  -0.662 0.508257    
onset_prev       -4.036e-04  5.529e-05  -7.300 4.01e-13 ***
slope_prev        1.152e-02  3.088e-03   3.732 0.000195 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0179 on 2178 degrees of freedom
  (184 observations deleted due to missingness)
Multiple R-squared:  0.6593,    Adjusted R-squared:  0.6567 
F-statistic:   248 on 17 and 2178 DF,  p-value: < 2.2e-16

CRF predicting mean difference

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_meandiff_socdist <- cforest(mean_diff_socdist ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n +
                                  airport_distance + republican + medage + male + popdens + 
                                  manufact + tourism + academics + medinc + physician_pc,
                           data = df_us_cpt_socdist,
                         controls = ctrls)

crf_meandiff_socdist_varimp <- varimp(crf_meandiff_socdist, nperm = 1)
crf_meandiff_socdist_varimp_cond <- varimp(crf_meandiff_socdist, conditional = T, nperm = 1)

crf_meandiff_socdist_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

crf_meandiff_socdist_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))